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ELASTIC  DEFORMATIONS  OF  A  ROTATING  SPHEROIDAL  EARTH 
DUE  TO  SURFACE  LOADS 


INTRODUCTION 

The  study  of  solid  Earth  dynamics  and  more  particularly  of  the  elastic  deformations 
of  the  Earth  has  been  in  the  foreground  of  investigation  for  a  couple  of  decades  and  has 
important  ramifications  into  geodesy,  geophysics,  and  astronomy.  Before  examining  some 
of  the  major  problems  related  to  an  elastic  Earth,  it  is  appropriate  to  furnish  a  few 
definitions  and  basic  facts. 

The  Earth’s  gravity  field  is  due  to  two  forces:  the  gravitational  attraction  exerted  by 
the  whole  Earth  at  each  point  according  to  Newton’s  law  and  the  centrifugal  force  due  to 
the  Earth’s  rotation.  The  magnitude  of  the  vector  field  is  the  intensity  of  gravity,  and  its 
direction  provides  the  direction  of  the  vertical  at  that  point. 

Next  in  order  of  importance  is  the  tidal  potential  generated  by  the  Moon  and  Sun 
(also  referred  to  as  the  luni-solar  potential).  This  potential  varies  with  time  because  of  the 
relative  motion  of  these  two  bodies  with  respect  to  the  Earth.  In  particular,  the  true  orbit 
of  the  Moon  is  quite  elaborate  and  differs  considerably  from  a  Keplerian  ellipse.  Various 
perturbing  terms  must  be  taken  into  account;  foremost  among  them  is  the  evection  term. 

The  most  important  consequence  of  these  tidal  perturbations  is  the  ocean  tide,  whereby 
the  free  surface  of  the  ocean  conforms  to  the  changing  field  and  remains  an  equipotential 
surface.  Another  consequence  of  this  tidal  potential  is  the  Earth’s  tide  or  the  elastic  defor¬ 
mation  of  the  Earth,  which  was  first  considered  by  Darwin  in  1883  [1]  and  ascertained  to 
be  about  a  third  of  the  ocean  tide. 

Our  knowledge  of  the  interior  of  the  Earth  is  provided  by  seismic  studies,  which 
studies  have  been  responsible  for  ascertaining  the  velocity  of  propagation  of  the  seismic 
waves  and  for  measuring  the  spheroidal  and  toroidal  oscillations  experienced  by  the  Earth 
subsequent  to  a  seismic  event.  From  this  knowledge  one  can  ultimately  determine  a 
rheological  profile  of  the  Earth  comprising  the  variation  of  its  density  and  of  the  Lame 
elastic  parameters  along  various  layers.  Love  introduced  in  1909  [2]  special  parameters 
which  bear  his  name  and  which  are  related  to  the  elastic  deformation  of  a  homogeneous 
spherical  Earth.  Like  the  moments  of  inertia,  these  Love  numbers  represent  intrinsic 
properties  of  the  planet  and  provide  estimates  for  the  maximum  rigidity  of  the  solid  inner 
core,  which  seismology  cannot  furnish,  since  no  transverse  wave  can  penetrate  the  liquid 
outer  core. 

Henceforth  in  this  report,  by  Earth  tides  we  shall  signify  the  elastic  deformations  of 
the  Earth  due  to  the  luni-solar  potential.  Let  us  now  examine  a  number  of  phenomena 
which  are  related  to  and  can  be  better  explained  through  the  concept  of  an  elastic  Earth: 
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•  Earth  tides  produce  a  periodic  variation  in  the  direction  of  the  vertical  which  is 
reflected  in  the  variation  of  latitude  for  an  astronomical  observatory  and  for  satellite 
tracking  stations.  Thus  Earth  tides  have  implications  to  observational  astronomy  and  the 
orbit  determination  for  artificial  satellites. 

•  Long-period  Earth  tides  perturb  the  Earth’s  axis  having  the  largest  moment  of 
inertia.  To  conserve  angular  momentum,  the  Earth  reacts  by  changing  its  rotational  rate. 
First  predicted  by  Jeffreys,  this  phenomenon  was  later  confirmed  by  the  use  of  atomic 
standards.  The  present  treatment  of  the  problem  is  still  based  on  the  assumption  of  a  rigid 
Earth  [3]. 


•  The  dissipation  of  energy  due  to  the  tidal  motion  contributes  to  the  secular 
retardation  of  the  Earth’s  rotation. 

•  The  luni-solar  tidal  force  acts  on  the  spheroidal  Earth  by  exerting  a  torque 
about  an  equatorial  axis.  The  rotating  Earth  reacts  as  a  gyroscope,  giving  rise  to  an  angular 
motion  which  can  be  decomposed  into  a  precession  of  the  Earth’s  axis  about  the  normal  to 
the  ecliptic  augmented  by  a  minor  oscillation  about  the  instantaneous  axis  known  as  a 
nutational  motion.  The  most  recent  treatment  of  this  phenomenon  is  due  to  Kinoshita  [4] 
and  also  assumes  a  rigid  Earth. 

•  The  presence  of  the  liquid  outer  core  within  the  Earth  produces  a  resonance  for 
oncoming  waves  having  a  tesseral  harmonic  character.  Poincare  [5]  first  demonstrated  that 
a  fluid  ellipsoid  contained  within  a  rotating  mantle  exhibits  a  free  nutational  mode.  The 
most  advanced  theory,  formulated  by  Molodenskii  [6] ,  ascertains  this  resonance,  but  it  still 
is  far  from  constituting  a  final  body  of  results  on  this  subject  (see  Melchior  [7] ). 

•  Amplitudes  of  the  Earth  tides  reach  as  much  as  0.3  to  0.4  m  and  cannot  be 
neglected  with  regard  to  lunar  and  satellite  laser  measurements,  which  have  accuracies  good 
to  few  centimeters. 

•  The  tidal  perturbations  of  the  Earth’s  potential  on  a  typical  satellite  orbit  can 
amount  to  50  m  and  must  be  taken  into  account. 

•  A  better  determination  of  the  surface  strains  produced  by  the  Earth’s  tides  will 
increase  the  accuracy  that  can  be  reached  in  long-base  interferometer  measurements  as  well 
as  by  other  instrumentation  which  depends  on  Earth-fixed  length  standards. 

•  Cubic  expansion  of  the  Earth’s  crust  will  produce  oscillations  in  water  wells 
and  can  have  a  bearing  in  ascertaining  the  storage  capacity  and  porosity  of  an  aquifer.  In 
general,  the  tidal  variation  of  underwater  fluids  can  play  a  factor  in  tectonic  processes  and 
in  their  diagnosis. 

Considering  the  variety  of  problems  where  the  influence  of  the  tidal  potential  plays  a 
significant  role,  it  is  appropriate  that  more  precise  and  sophisticated  methods  be  developed 
to  ascertain  both  the  ocean  and  Earth  tides.  The  present  knowledge  of  the  global  ocean 
tides  is  based  on  the  numerical  integration  of  the  Laplace  tidal  equations.  However,  the 
simplifying  assumptions  of  a  rigid  Earth  and  of  a  homogeneous  incompressible  ocean  in 
hydrostatic  equilibrium  which  are  used  to  solve  those  equations  seem  to  vitiate  the  results 
and  fail  to  give  close  agreement  with  tidal  observations  at  midocean  islands.  The  problem 
should  be  formulated  as  a  boundary-value  problem  by  considering  the  continental  boundary 
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of  the  ocean  basin  and  considering  the  ocean  bottom  with  friction  losses.  Furthermore,  the 
attraction  of  the  ocean  mass  should  also  be  taken  into  account,  and  this  will  transform  the 
Laplace  tidal  equations  into  a  system  of  integro-differential  equations. 

The  Earth  tides  are  due  primarily  to  two  effects.  One  is  the  yielding  of  the  elastic 
Earth  to  the  gravitational  attraction  of  the  Moon  and  the  Sun.  This  direct  effect  goes  by 
the  name  of  body  tide  and  is  rather  well  understood,  because  astronomical  observations 
have  furnished  precise  information  on  the  forcing  function.  There  is  also  an  indirect  effect 
due  to  the  yielding  of  the  Earth  to  the  load  of  the  harmonically  varying  ocean  tide  upon  the 
continents;  this  is  also  known  as  load  tide.  The  corresponding  forcing  function  is,  however, 
not  well  known,  not  only  because  of  our  scanty  information  of  the  ocean  tides  but  also 
because  we  still  do  not  know  the  influence  played  by  the  structural  heterogeneities  in  the 
Earth’s  crust  and  upper  mantle.  However,  it  is  evident  that  load  tides  and  ocean  tides  are 
interrelated  phenomena. 

The  direct  and  indirect  Earth  tides  are  both  ultimately  derived  from  the  same  astro¬ 
nomical  input;  consequently  they  will  have  the  same  frequencies.  An  eventual  lag  of  these 
tides  with  respect  to  the  acting  potential  should  give  information  about  the  viscosity  of  the 
Earth.  However,  the  spatial  behavior  of  these  two  tides  is  quite  different :  the  body  tide 
varies  smoothly  over  the  Earth,  whereas  the  load  tide  presents  irregularities  because  of  the 
discontinuity  in  the  forcing  function  at  the  coastline  and  also  because  the  ocean  tide  exhibits 
a  circulation  motion  around  the  amphidromes. 

Crustal  loading  constitutes  a  phenomenon  whose  frequencies  are  intermediate  between 
those  of  seismic-wave  transmission,  which  is  an  elastic  event,  and  tectonic  deformations, 
which  have  an  anelastic  character.  Thus  one  might  expect  that  the  study  of  loading  defor¬ 
mations  are  useful  in  determining  the  threshold  of  anelasticity  for  the  crust  and  upper 
mantle  layers. 

Ocean  load  effects  are  responsible  in  part  for  the  variation  of  the  vertical  component 
of  the  gravity  field  (also  referred  to  as  tidal  gravity),  for  surface  strains  (strain  tide),  and  for 
the  deflection  of  the  vertical  (tilt  tide).  As  a  consequence,  if  for  a  certain  region  the  Earth 
model  obtainable  from  seismic  evidence  is  better  known  than  the  cotidal  charts  of  the  ocean 
tide,  then  one  could  expect  [8]  that  accurate  tidal  gravity  measurements  on  adjacent  lands 
(which  now  can  reach  an  accuracy  of  1  microgal)  can  be  used  to  improve  our  knowledge  of 
the  ocean  tide.  Conversely,  should  the  reverse  be  true,  the  ocean  tide  can  provide  some 
information  on  the  elastic  properties  of  the  crust. 

The  purpose  of  this  report  and  possibly  of  future  reports  and  papers  is  to  provide  a 
comprehensive  study  of  the  load  tides  for  a  rotating  spheroidal  Earth.  The  basic  work 
which  constitutes  our  theoretical  background  goes  back  to  Longman  [9,10]  and 
Farrell  [  11  ] .  These  two  authors  have  provided  the  basic  equations  and  boundary  conditions 
for  the  load  tide  of  an  elastic  nonrotating  spherical  Earth.  More  recently,  Smith  [12]  ob¬ 
tained  a  model  for  the  elastic  deformations  of  a  rotating  spheroidal  Earth.  His  formulation 
is  very  general  and  mathematically  correct;  however,  it  was  found  to  be  rather  cumbersome 
to  attack  from  a  computational  point  of  view  due  to  his  use  of  a  special  set  of  harmonic 
functions. 

For  future  geophysical  developments,  the  effect  of  the  Earth’s  rotation  on  the  Earth 
tide  should  be  studied  and  ascertained  from  a  numerical  point  of  view.  We  have  therefore 
approached  here  the  same  problem  of  determining  the  elastic  deformations  of  a  rotating 
spheroidal  Earth  and  have  been  able  to  provide  a  much  simpler  analytic  representation  of 
such  deformations.  Our  success  is  primarily  due  to  our  having  been  able  to  write  the 
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equations  pertaining  to  the  nth-order  harmonic  based  on  a  simple  decomposition  of  the 
product  of  Legendre  polynomials  and  their  derivatives  into  spheroidal  components.  This 
simple  representation  of  the  governing  equations  augurs  well  for  their  numerical  solution. 

We  present  in  this  report  the  theoretical  work  that  leads  to  the  representation  of  the  equa¬ 
tions  valid  for  the  deformation  of  a  rotating  spheroidal  Earth  and  also  give  an  insight  on 
future  numerical  work.  The  full  expansion  of  this  numerical  endeavor  will,  however,  consti¬ 
tute  material  for  a  future  report. 


LINEARIZED  NAVIER-STOKES  EQUATIONS 


With  respect  to  an  inertial  frame  (O;  xltx2,  x3),  with  the  origin  O  coinciding  with  the 
centroid  of  the  planet,  the  equations  representing  the  displacement  field  u  from  the  initial 
reference  state  can  be  written 


d2u(  _  ay  bo‘i 
P  dt2  P  d*;  +  dxy  ’ 


(1) 


where  the  summation  convention  applies  to  repeated  indices.  Here  the  quantities  u,-  are  the 
components  of  the  displacement  field  u  and  the  quantity  p  is  the  mass  density  and  is  the 
sum  of  pQ,  the  mass  density  of  the  reference  state,  and  py,  the  change  in  mass  density  due 
to  the  displacement  field,  with  py  and  the  components  u,  being  related  by  the  continuity 
equation 


Pi  =  -  ^7  (Po«l>  •  <2> 

Further  in  Eq.  (1),  V  is  the  total  potential  and  is  the  sum  of  V0,  the  potential  of  the  refer¬ 
ence  state,  and  Vr1,  the  potential  due  to  the  deformation.  V0  consists  of  the  gravitational 
potential  and  the  rotational  potential;  it  satisfies  the  Poisson  equation 


V2V0  =  -4irGp0  +  2u>2  ,  (3) 

where  G  i^the  gravitational  constant  and  w  is  the  Earth’s  angular  velocity.  Vy  satisfies  also 
a  Poisson  equation 


V2^  =  -4ir Gpy  . 


The  components  of  gravity  are  then 


&0,i 


9Vo 

9x,- 


8\,i 


dxj 


The  gravitational  attraction  is  given  by 


,  ,  4irG  r  .  .  2  j  4irG  —  .  , 

v>  =  J  po(r)r  dr  =  “  rp0(r)  , 


(4) 


(5) 


(6) 
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where  r  is  the  radial  distance  measured  from  the  centroid  and 


P0(r)  =  -j  f  P0(r)r2  dr 
r 


is  the  mean  density.  The  derivative  of  the  gravitational  attraction  with  respect  to  r  satisfies 
the  relation 

7o(r)  =  -  p0  +  4t rGp0  .  (8 

Finally  in  Eq.  (1),  a,,  are  the  components  of  the  total  stress  field.  This  total  field  consists 
of  two  fields.  The  first  is  the  initial  stress  field,  which  we  assume  to  be  hydrostatic  equilib¬ 
rium  and  whose  components  T ):  satisfy  the  conditions 


uy0 

T ij  ~  -Po^ij  >  ~  ’ 


where  pQ  is  the  hydrostatic  pressure  and  the  quantities  6 are  the  Kronecker  deltas.  The 
second  is  the  stress  field  valid  for  a  perfectly  elastic  and  isotropic  medium,  whose  compo¬ 
nents  Tjj  are  related  to  the  Lame  elastic  parameters  X  and  p  by  the  relation 

3  uk  /du{  <*Uj\ 


Since  we  are  assuming  infinitesimal  elastic  deformations,  the  above  relationship  applies  to 
the  coordinates  that  the  points  of  the  medium  had  before  the  occurrence  of  the 
deformation.  The  same  assumption,  however,  cannot  be  made  for  the  initial  stress  field. 

As  a  consequence,  the  total  stress  at  the  undeformed  points  of  the  medium  must  be 
written  as 

ZTij 

°ij  T‘j  dxk  Uk  +  T‘j  '  (11 

By  differentiating  Eq.  (9),  we  can  easily  get 

^  Tij  dp0 

dx k  =  ~dx^  S|>  =  _p0^0,fe6,y>  (li 

and  this  can  be  used  to  rewrite  the  fundamental  equations  of  the  displacement  (Eq.  (1))  as 


Po*l,i  +  PlSo,i  +  foT  ( PoSo,kuk )  + 


where  we  have  neglected  the  product  p^gx  ,  because  it  is  of  the  second  order  in  the 
displacements. 
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In  a  previous  work  [13] ,  we  gave  details  on  how  to  write  the  previous  equation  in 
vectorial  form,  and  we  obtained 

P  =  Po?i  +  Pi?o  +  V(p0u-I0)  +  E(\,u),  (14) 

dt- 

where 

E(X,  p)  =  (X  +  2  p)V(V-  u)-pVxVxu  +  (VX  +  Vp)(V  •  a) 

+V(u  •  Vp)  +  V  X  (u  X  Vp)  -  uV2p  (15) 


stands  for  the  agglomerate  of  terms  due  to  the  elastic  parameters  X  and  p. 


The  first  three  terms  on  the  right-hand  side  of  Eq.  (14)  represent  respectively  the 
action  of  the  perturbed  gravity  field  on  the  initial  density  profile,  the  action  of  the  initial 
gravity  field  on  the  deformed  density,  and  the  contribution  due  to  the  work  performed  by 
the  elastic  displacement  vector  against  the  initial  hydrostatic  field.  Equation  (14)  is  essen¬ 
tially  a  linearized  version  of  the  Navier-Stokes  equation  when  allowance  is  made  for  the 
variation  of  the  material’s  elastic  parameters  and  quadratic  terms  in  the  displacements 
are  dropped. 


The  acceleration  of  u  with  respect  to  the  inertial  frame  is  given  by 

d2u  d2u  du)  ->  -►  2t* 

- =  — —  +  2co  X  —  +  -77-  X  u  +  (to  •  u)u>  -  c o^u  , 

dt 2  dt2  dt 


(16) 


where  the  partial-derivative  symbol  refers  to  the  variation  of  a  vector  with  respect  to  the 
rotating  frame. 

In  this  report,  we  are  interested  in  ascertaining  permanent  deformations  of  a  rotating 
Earth,  so  that  the  first  and  second  time  derivatives  of  u  with  respect  to  the  rotating  frame 
must  vanish.  Also,  we  are  assuming  a  constant  angular  rotation  and  no  oscillation  for  the 
position  of  the  rotational  axis,  which  is  equivalent  to  saying  that  the  time  derivative  of  c5 
is  also  vanishing.  We  are  therefore  left  with  the  last  two  terms  of  Eq.  (16),  which  pertain 
to  the  centrifugal  acceleration. 

Reverting  to  Eq.  (14),  we  shall  write 

P0[(c2  •  u) c3  -  u>2u]  =  p0gx  +  Pig0  +  V(p0u-  g0)  +  p)  ,  (17) 

where  the  term  p^u  has  been  neglected,  being  of  the  second  order  in  the  displacements. 

This  is  the  fundamental  equation  of  motion,  and  we  shall  seek  solutions  to  this  equation  in 
the  form  of  spheroidal  and  toroidal  deformations. 

In  the  sequel  we  shall  assume  rotational  symmetry,  whereby  the  displacement  vector  u 
depends  only  on  the  radial  distance  r  and  the  colatitude  0  of  the  position.  The  spherical 
harmonics  Yn(9,  0)  will  then  reduce  to  the  Legendre  polynomials  P„(cos  9).  We  use  spher¬ 
ical  coordinates  (r,  9,  0)  with  the  origin  at  the  centroid,  9  being  the  colatitude  and  0  the 
longitude.  Accordingly,  the  rotational  vector  is  representable  as 
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OJ  s  {cocos  0  ;  -co  sin  0  ;  0}  (18) 

and  the  unperturbed  gravity  vector  as 

[av0  i  dvQ  } 

io  »  w0(r,  8).  j-jij  oj  .  (19) 

SPHEROIDAL  DEFORMATIONS 

The  spheroidal  displacement  vector  u  has  components  of  the  form 


u  ~  2]  U„(r)Pn{c°s  6) , 
n*0 

u  =  Z  K(r)  Tf .  (20) 

n  =  l 
UJ  =  0 

and  has  the  property  that  the  radial  component  of  its  curl  vanishes.  We  shall  determine  the 
two  sequences  of  unknown  functions  Un(r )  and  Vn(r)  by  imposing  the  conditions  that  the 
two  series  in  Eq.  (20)  satisfy  the  fundamental  equations  of  motion:  Eq.  (17). 

Another  quantity  of  physical  significance  related  to  our  problem  is  the  deformed 
potential  which  we  assume  to  be  of  the  form 


Vj(r,  0)  =  ]2  Rn(r)Pn  (cos  0)  .  (21) 

ii=0 

This  introduces  another  sequence  Rn(r)  of  unknown  functions.  It  follows  that  g1, 
being  the  gradient  of  is  given  by 


§1  =  Wj 


RnPn 


R n  K x 
r  dd 


Henceforth  primes  shall  denote  derivatives  with  respect  to  the  radial  distance  r. 


(22) 


From  Eq.  (20),  we  obtain  that  the  dilatation  A  can  be  written  as 


A  =  V  •  u  =  Xn(r)Pn(cosQ)  , 
n  =  0 


(23) 
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( 

1 


where 

Xn  =  V'n  +7Un~  ~T  1 

and  that  the  perturbed  density,  as  given  by  the  continuity  equation,  is 

oo 

Pi  =  ~V‘(P0u)  =  (Unp'0  +  X„p0)P„(cos  d)  . 

n=0 

The  last  of  the  inertia  terms  appearing  in  Eq.  (17)  can  be  developed  as 

V(p0«-fb>  =  *?o)VPo  +  PqV(»  •  g0) . 


(24 


(25. 


(26) 


where  p0  is  a  given  function  of  r.  The  above  expression  implies  the  knowledge  of  the  partial 
derivatives  of  V’0  up  to  and  including  the  second  order.  V0  refers  to  the  initial  hydrostatic 
regime.  In  Appendix  A  we  have  expressed  these  derivatives  in  terms  of  the  rotational 
deformation  obtainable  from  the  Clairaut  equation. 

To  complete  the  analytical  derivation  of  the  equations  of  motion,  one  first  must 
eliminate  the  trigonometric  terms  via  the  Legendre  polynomials.  The  required  identities 
simply  reduce  to 


sin2  9  =  |-  (1  -P2), 
cos2  0  =  ~  (1  +2P2)  , 


sin  0  cos  0 


I  dA 

3  bO  ' 


(27) 


Then  one  must  express  the  spheroidal  components  of  the  products  P2^t»  dPn/dO, 
PndP2  /30,  and  (dP2/dO)  •  (dPn/dff)  for  a  generic  value  of  n.  This  is  equivalent  to  writing 
these  products  as  linear  combinations  of  Pn  and  9P„/00. 


This  operation  has  been  accomplished  in  Appendix  B  by  elaborating  on  certain  well- 
known  relationships  valid  for  the  Legendre  polynomials.  Let  us  also  note  that  the  second 
derivatives  of  the  P’s  can  be  eliminated  by  means  of  the  fundamental  equation  of  definition 

^sin  0  =  ~n(n  +  1)  sin  0  Pn  .  (28) 

The  elastic  terms  E( A,  p)  which  are  given  by  Eq.  (15)  have  also  been  expanded  and  found  to 
agree  with  the  results  of  Alterman  et  al.  { 14]  and  of  Longman  [9] .  We  therefore  shall 
refrain  from  discussing  their  derivation  in  this  report. 
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The  radial  component  of  Eq.  (17)  was  found  to  be  as  follows: 

y0p0(4  -  3D0)Un  -  n(n  +  l)y0pQVn  +  rp0Rn  +  r(\Xn  +  2pU'n )' 

+  7  WK  -  4Un  +  n(n  +  l)(3Vn  -  rV ;  -  t/„)l 

+  rp0GJ2(A1Un-2  +  A2  Un  +  A3Un  +  2  +  A4Vn-2  +  A5  Vn  +  ^6^0  +  2 

+A7rVn-2  +  A8rVn  +  A9rVn*2)  =  (29) 

The  two  transversal  components  of  the  same  equation  give  identical  contributions,  which 
amount  to 


-70P0 Un  +  PO Rn  +  +  r  [  M  [k  ~  ~  +  ^ 

+  7  [5£/„  -  -  2n(n  +  l)y„  +  3rV'] 

+  rp0W 2  (-81^-2  +  ^2^n  +-®3^n  +  2  +^4^n-2 

+^5^n  +B6Vn  +  2  +  B7rUn-2  +  B8rUn  +  B9rUn  +  2)  =  0  •  (30) 

Here  70(r)  is  the  gravitational  attraction  given  by  Eq.  (6)  and  D0  =  p0/p0  is  the  ratio 
between  density  and  mean  density  as  defined  by  Eq.  (7).  The  quantities  An  and  Bn 
(n  =  1,  2,  ....  9)  are  two  sets  of  known  functions  which  depend  on  the  radial  distance  r, 
the  density  p0,  the  hydrostatic  deformation  f2i,  and  their  derivatives  p'0  and  f2x.  They 
are  related  to  the  quantities  A*,  B*,  and  C*  obtained  in  Appendix  A.  It  is  premature  at 
this  stage  to  explicitly  write  down  these  18  A  and  B  functions,  simply  because  for  purposes 
of  numerical  evaluation  it  is  more  expedient  to  deal  with  particular  linear  combinations  of 
the  A  and  B  functions,  and  these  will  be  furnished  later  by  Tables  1  and  2. 


We  need  an  extra  equation  to  complete  the  determination  of  the  three  sets  Un,  Vn, 
and  Rn  of  unknown  functions;  this  is  provided  by  the  Poisson  equation  (Eq.  (4)),  which 
explicitly  becomes 

r2R"n  +  2 rRn  -  n(n  +  l)Rn  =  4nGr2(p'0U„  +  p0Xn)  .  (31) 

The  terms  appearing  within  Eqs.  (29)  and  (30)  and  which  are  due  to  the  rotation  to  relate 
to  the  coefficients  of  the  harmonics  of  order  n  -  2,  n,  and  n  +  2.  Thus  we  have  an  infinite 
set  of  equations  relating  three  harmonics. 

Related  to  these  three  sets  of  unknowns  are  the  Love  numbers,  which  are  nondimen- 
sional  quantities  defined  as 
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1  +  *„(»■)  =  Rn(r)lnQ(.r)  . 

hn(r)  =  Un(r)/r,  (32) 

*n(r)  =  V„(r)/r  . 

They  depend  on  the  order  n  and  the  radial  distance  r. 

REDUCTION  TO  NORMAL  FORM 

The  numerical  evaluation  of  the  second-order  differential  Eqs.  (29),  (30),  and  (31) 
requires  that  they  be  reduced  to  their  normal  form,  that  is,  to  the  first-order  system 
y'  =  F(x,  y),  which  is  amenable  to  programming  on  a  computer.  This  reduction  can  be 
achieved  by  introducing  in  addition  to  the  original  set  of  variables 

Un(r)  (radial  displacement), 

Vn(r)  (transversal  displacement),  (33) 

Rn(r )  (change  in  gravitational  potential) 
a  new  set  of  three  variables.  The  first  new  variable  is 
En(r)  =  XAn  +  2pf/„' 

=  (X  +  2p)U„  Un  -  /i(«+l)£  V„,  (34) 

which  appears  in  the  expansion  of  the  radial  stress 

OO 

rrr  =  2^  En{r)Pn(cos  0)  .  (35) 

n  =  0 

The  second  new  variable  is 

W  =  »(<  ~TVn  +TUn )-  (W 

which  is  connected  with  the  transversal  stress 

wLww-  (37) 

n~l 

The  third  new  variable  is 

Hn(r)  =  R'n  -  4t rGp0Un  ,  (38) 

which  represents  the  change  in  gravitational  flux. 
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The  linearized  Navier-Stokes  equations  representing  the  spheroidal  deformations  can 
now  be  written  down  as  the  following  differential  system  of  six  equations: 

U*  =  "  X  +  2(i  7  Un  +  n(n  +  1)  X  +  2fi  ~  Vn  +  X  +  2p  En  ’ 

<  -  -  7  ^  +  7  +  . 

fj;  =  47rGp0t/„  +  Hn  , 


r,i  n(n  + 1)  n(n  +  1)  2 

*»  =  _47rG^o  — -  ^  + - -  Rn~yHn, 


K  =  4 


ToPo  /  3X  +  2/A  /Lt 


+n(n  +  1) 


X  +  2p  )  r2  1  U" 
7qPq  /  3X  +  2(i\  2p 


X  +  2 At  /  r2 


4(i  i  „  n(n+l) 

T  E„  + - ^  F„  -  p^ 


X  +  2(i  r  « 

+p0co2fC1i/„.2  +  C2f/„  +  C3Un  +  2  +  C4F„_2  +  c5vn  +  C6V„+2 

+  7  (G7^n-2+G8^+^  +  2)], 

p,  f 70Po  /  3X  +  2n\2n 
'"'L  r 


2(i 

—  V 
r2  n 


(39) 


X  +  2p  )  r2  '  U* 

(2  n2  +  2n~  1)X  +  2(n2  +  n-  1  )(i 
X  +  2(i 

_  X  1  _  3 

r  "  X  +  2(i  r  "  r 

+p0w2[D1Un_2  +I>2l/n  +D3Un  +  2  +D4vn.2  +  DhVn  +  Z>6Vn  +  2 
+  X_+_2(l  ^7^1-2  +  D%En  +  •09^+2)!  • 

Here  primes  denote  as  usual  derivatives  with  respect  to  the  radial  distance  r.  The  two  sets 
of  dimensionless  functions  C  and  D  are  explicitly  given  in  Tables  1  and  2  respectively  in 
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terms  of  the  A*  B *,  and  C*  quantities  that  were  obtained  in  Appendix  A.  They  ultimately 
depend  on  the  solutions  of  the  Clairaut  equation. 

Before  proceeding  to  the  numerical  solutions  of  Eqs.  (39),  we  must: 

•  Determine  the  conditions  that  must  be  satisfied  in  order  that  these  equations 
have  regular  solutions  at  the  center  of  the  configuration  (r  =  0),  that  is,  solutions  which  are 
expressible  as  infinite  power  series  about  the  center; 

•  Ascertain  the  limiting  form  of  these  equations  for  the  case  n  =  0  (zero-order 
harmonic),  since  the  nth-order  harmonic  depends  on  the  harmonic  of  order  n  -  2; 

•  Ascertain  the  limiting  form  of  these  equations  which  is  valid  within  the  liquid 
outer  core  where  p  =  0 ;  and 

•  Determine  the  conditions  that  the  system  must  satisfy  at  the  free  surface 
(r  =  a j )  depending  on  the  loading  conditions. 


REGULARITY  CONDITIONS  FOR  THE  SOLUTIONS  AT  THE  CENTER 
WHEN  ju  *  0  AND  n^O 


We  assume  power -series  expansions  in  the  neighborhood  of  r  =  0  for  each  of  the 
six  variables 


Un(r)  =  Un0  +  Unlr  +  C7n2r2  + 
Hn(r)  =  Hn0  +  Hnl  r  +  Hn2r  2  +  - 


(40) 


and  replace  them  within  both  sides  of  Eqs.  (39).  From  a  physical  point  of  view,  we  must 
also  require  that  the  components  of  the  displacement  be  zero  at  the  origin : 

^0=0,  V„o  =  0  •  <«) 

Since  the  right-hand  sides  of  the  equations  contain  terms  in  l/r2  and  1/r,  the  regularity 
conditions  will  be  obtained  by  imposing  the  vanishing  of  the  coefficients  of  the  two 
negative  powers  of  r  and  by  equating  the  constant  terms  appearing  on  both  sides.  All  other 
positive  powers  of  r  do  not  yield  any  extra  conditions,  since  we  deal  with  the  limit  when  r 
approaches  zero. 

Let  us  remark  that  from  Eq.  (6)  it  follows  that 


lim 

r-0 


y0(r)/r  =  p0(0) 


(42) 


is  finite  and  also  that  the  rotational  terms  will  not  give  any  contribution,  since  they  are 
formed  with  the  U  and  V  terms. 
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Each  equation  can  give  rise  to  at  most  three  conditions.  We  ultimately  arrive  at  ten 
conditions  among  the  coefficients  of  the  variables.  Five  of  these  conditions  can  be 
immediately  analyzed: 


^n0 

=  o, 

Rnl 

ii 

o 

(n  -  l)(n  +  2)Hn0 

=  0, 

(43) 

Uni 

-  i  F 
"  M  n0  ’ 

Rn2  -  n(n  +1)  Hnl  +  47rGPoKil  ' 


From  the  second  and  third  conditions,  we  see  that  when  n  1 ,  we  must  have 

Hn0  =  Rnl  =  0,  n  *  1,  (44) 

whereas  when  n  =  1 ,  one  can  only  say 


*H  =  H10  .  (45) 

Next,  we  have  three  linear  homogeneous  relations  among  ^n0’  ^nO’  and  The  determi¬ 
nant  made  with  the  coefficients  of  the  unknowns  vanishes  for  n  =  2.  Thus  in  this  case  we 
reach  solutions  of  the  form 


V21  =  ^20  >  E20  =  2F20  , 

whereas  when  n  =£  2,  there  is  only  the  solution 

^nl  =  En0  =  En0  =  0  > 


(46) 


(47) 


and  this  leads  to  t/nl  =  0. 


The  last  two  conditions  of  the  set  relate  the  six  quantities  Un2,  Vn2 ,  Enl,  Fnl,  Rnl, 
Hn0.  In  the  case  ni=  1,  when  Hn0  =  Rnl  =  0,  these  conditions  can  be  used  to  ascertain 
Un2  and  Vn2  in  terms  of  En\  and  Fn-y  according  to  the  formulas 


2ju(n  -  l)(n  +  2)  Vn2 


(48) 


2ji(3X  +  2 y.)(n  -  l)(n  +  2 )Un2  =  [(5n2  +  5rt  -  1)X  +  6(n2  +  n  -  l)q]Efll 


-n(n  +  l)[(2n2  +  2n  -  13)X  +  2(n2  +n-  5)^]F„i . 
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When  n  =  1 ,  the  same  conditions  can  be  used  to  ascertain  E , ,  and  F,  i  in  terms  of  t/,  „  V,  0 
and  Hl0  =  according  to  11  12>  12’ 


„  „  /  3X  +  2/i\  X  +  2u 

£l1  ' 2,1  (rrjr)  (t,«  - v«>  -  rrv  «o»io  • 

/3X  +  2(j\  u 

11  "  ITT  4p  j  (f/l2  "  Vl2)  '  r+  4p  PoHlO  ■ 


Consequently,  for  n  -  1  the  series  expansions  are 


(49) 


Ux  = 


*1  = 
= 


**12 ' 
^i2  r2 


+ 

+ 


H10r  +  (3/2)f/ur2  + 


(50) 


Fur  + 


Fur  + 


#1  =  **10  +  **llr  +  '"> 

with  H10,  //n,  f/12,  and  V\2  arbitrarily  chosen  and  with  Fn  and  Fn  given  by  Eqs.  (49). 
When  n  =  2,  the  series  expansions  are 


**2  = 
^2  " 


*2  = 


(*20  /M)r  +  ^22r2 
(F20/2jU)r  +  V22r2 


+ 

+ 


(1/2)  (tf21  +  F20 


■®2  2F20  +  F21r 

*2  =  *20  +  *21 r 


r*  + 


**2  = 


**21 ' 


+ 

+ 

+ 


(51) 


with  arbitrary  values  for  F20,  F21,  F21,  and  //21  and  with  U22  and  V22  given  by  Eqs.  (48). 
Finally,  when  n>  3,  we  can  write  expansions  of  the  sort 


Un  = 


K  = 


Rn  = 


**„2r2 

v^2 


+ 

+ 


[3/n(n  +  l)]ffnlr2  + 


(52) 


En\r  + 


Fn  =  *nl r  +  — , 
**n  =  **„ir  +  -, 
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with  Enl,  Fni,  and  Hnl  arbitrarily  chosen  and  with  Un2  and  Vn2  also  determined  by 
Eqs.  (48). 


FUNDAMENTAL  EQUATIONS  OF  MOTION  FOR  THE  ZERO-ORDER  HARMONIC 

When  n  =  0,  there  is  no  transverse  spheroidal  component  for  the  displacement  and  for 
the  stress,  since  dP0lbd  =  0.  This  translates  into  the  fact  that  there  exists  no  equation  for 
the  6  component:  U0(r)  is  arbitrary  and  so  is  F0(r),  which  depends  on  V0(r).  The  system 
of  equations  in  its  normal  form  cannot  contain  any  equation  pertaining  to  V0  and  F0;  these 
two  variables  do  not  appear  within  the  remaining  four  equations  relating  UQ,  Rq,  Eq,  and 
H0.  From 


Ho  =  -<2/r)H0 


(53) 


we  solve  and  get  H0(r)  =  h/r2.  Regularity  at  r  =  0,  however,  requires  that  h  must  be  taken 
equal  to  zero;  thus  H0  =  0.  The  remaining  equations  are 

rj'  _  _  2X  —  u  +  ^  p 

U°  X  +  2  n  r  U0  X  +  2/1  E°  ’ 

R'o  =  4tt Gp0U0  ,  (54) 


E 


0 


/3X  +  2/i  \  p 
\  \  +  2  fi)  T2  . 


UQ  - 


4/i 

X  +  2  p 


+p0cj2  [ C3(0) U2  +  C6(0) V2  +~C9(0)F2}. 
Regularity  at  r  =  0  imposes  the  conditions 

^00  =  *oi  =  ® 


(55) 


and 


E00  =  (3X  +  2p)Uqi  , 

/  3X  +  2  p\ 

E°'  =  (  X  +  6 p)  U°2' 

Thus  the  series  expansions  in  the  neighborhood  of  r  =  0  can  be  written  as 

U0  =  U01r  +  U02r 2  + 

*0  =  *00  +  -*02r2  +  ■- 
*0  =  *00  +  *01 r  +  ■"> 


(56) 


(57) 
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where  t/01,  UQ2,  R00,  and  RQ2  can  be  chosen  arbitrarily  and  E00  and  £01  are  obtainable 
from  Eqs,  (56). 

Note  that  V0(r)  anc*  F0(r)  are  arbitrary  functions.  They  do  not  appear  in  Eq.  (54) 
pertaining  to  n  =  0;  neither  will  they  appear  in  the  other  possible  case,  that  is,  among  the 
rotational  terms  of  Eqs.  (39)  when  n  =  2.  This  is  so  because  one  can  verify  from  Tables  1 
and  2  that  C4(2)  =  C7(2)  =  Z)4(2)  =  0. 


EQUATIONS  VALID  FOR  THE  OUTER  CORE 

To  investigate  the  equations  of  motion  valid  within  the  liquid  outer  core,  one  must 
evaluate  the  limit  reached  by  the  system  of  Eqs.  (39)  when  the  rigidity  p  approaches  zero. 

The  U  equation  becomes  simply 


U„ 


TUn  + 


n(n  +  1) 


V  +  —  E 
n  X  " 


The  V  equation,  when  written  as 

rFn  =  lim  p(rV„'  +  Un-Vn), 

P-0 


(58) 


(59) 


reveals  that  for  the  values  of  r  +  0  it  must  be  Fn  =  0,  whereas  the  terms  appearing  within 
the  parenthesis  could  be  arbitrarily  chosen. 

The  R  and  H  equations  do  not  contain  the  parameter  p  and  do  not  have  to  be  modified. 
The  E  equation  yields 


E 


n 


+  n  (n  +  1) 


TqPq 

r 


Vn  -  PoHn 


+p0co2[Cit4].2  +  C2Un  +  C3Un  +  2  +  C4V„_2  +  C5V,, 

+  C6V*  +  2  +  C7(rV^.2  +  Un_ 2  -  Vn_2) 

+C8(r^  +  Un  -  Vn)  +  C9(rV'n  +  2  +  Un  +  2  -  V„  +  2)l  •  (60) 

Here  use  has  been  made  of  Eq.  (59)  to  eliminate  the  terms  (r/p)Fp  in  favor  of  the  expression 
rV„  +  Up  -  Vp  (where  p  takes  on  the  three  values  n  -  2,  n,  and  n  f  2),  which  is  not  an 
indeterminate  form. 

The  F  equation  degenerates  into  an  algebraic  equation  and  can  be  used  as  a  checkout 
condition: 


En  =  yoPoU*  -  PoRn  +  rPo“2  [^>1^-2  +  °2 Un 
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Equation  (60)  is  a  differential  equation  not  in  its  normal  form,  since  both  En  and 
p0co2V'(p  =  n  -  2,  n,  n  +  2)  appear  in  it.  One  could  use  Eqs.  (60)  and  (61)  for  a  simulta¬ 
neous  determination  of  En  and  Vn.  Computationally  this  procedure  will  not  be  convenient, 
since  it  is  not  straightforward.  It  will  be  numerically  more  expedient  and  still  physically 
plausible  to  assume  Vn  =  0  throughout  the  whole  interior  of  the  outer  liquid  core  (to 
assume  the  vanishing  of  the  tangential  displacement),  since  Vn  has  been  left  arbitrary  in  the 
previous  considerations.  If  so,  and  for  the  case  n  ^  0,  we  get  the  following  system: 

Vn  =0-  Fn  =  0. 

K  =  -f  Un 

R'n  =  4 nGp0U„  +  Hn  , 

(62) 

,  n(n  +  1)  2 

=-V~  *n 


En  =  -4 


Un  -  P0Hn  +  P0w2t(Cx  +  C7)t/„_2 


+  (C2  +  Cs)Un  +  (C3+C9)[/n  +  2]  , 

=  yoPoUn  -  PoRn  +  [D\ Un-2  +  ^2^  +  D3Un  +  2 

+  ^  2  +  ^SEn  +  D9En  +  2^  • 

No  series  expansions  are  required  for  these  equations,  because  they  are  valid  within  the 
outer  core,  and  their  initial  conditions  at  a  value  of  r  0  are  furnished  by  the  final  values  of 
the  variables  at  the  inner  core  interface. 

Let  us  now  consider  the  more  particular  case  when  both  p  and  n  are  zero.  We  should 
then  ascertain  the  limit  of  Eqs.  (54)  when  the  rigidity  approaches  zero  and  eliminate  the 
term  (r/p)F2  according  to  Eq.  (59).  We  should  also  take  V2  =  0  in  agreement  with  Eq.  (62). 
The  two  functions  F0  and  V0,  which  in  Eqs.  (54)  were  left  arbitrary,  can  now  be  chosen  as 
both  vanishing.  We  are  left  then  with  the  system 

Uq  =  U0  +  l  E0  , 

Rq  =  AkGpqUq  , 


£o  =  -4 


U0  +  p0oj2[C3(0)  +  C9(0)]{/2  , 


V0  =  F0  =  H0  s  0  . 


t  - 
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If  we  differentiate  the  last  equation  of  Eqs.  (62)  with  respect  to  r  and  combine  it  with  the 
equation  appearing  prior  to  it  in  the  same  set,  we  get,  after  performing  a  number  of  simpli¬ 
fications,  a  relationship  which  symbolically  can  be  abbreviated  as 


{ypPo 


Pp 

PQ 


En  = 


P  0^ 


L 


r2  , 

4-  E 


p  =  n  -  2,  n,  n  +  2.  (64) 


The  vanishing  of  the  term 

7qPq  Po 
X  +  Po 


(65) 


has  been  known  as  the  Adams-Williamson  condition  [15] .  It  was  first  connected  to  this 
problem  by  Longman  [10]  and  was  afterward  extensively  studied  by  Pekeris  and 
Accad  [16].  In  the  absence  of  rotational  terms,  two  possibilities  were  available:  either 
the  Adams-Williamson  condition  should  have  been  valid  within  the  outer  core,  and  this  was 
found  to  be  unsatisfactory  from  a  physical  point  of  view,  or  one  should  have  En  =  0,  which 
according  to  Eqs.  (34)  and  (23)  amounts  to  zero  dilatation.  This  latter  condition  should  be 
discarded,  because  it  restricts  the  number  of  parameters  available  to  match  the  boundary 
conditions.  Thus,  the  inclusion  of  rotational  terms  in  the  study  of  elastic  deformation  of 
the  Earth  not  only  gives  rise  to  a  more  realistic  physical  model  but  also  avoids  the  use  of 
the  Adams-Williamson  condition,  which  use  has  been  controversial. 

Let  us  finally  note  that  the  Brunt-Vaisala  frequency  A  for  a  stratified  liquid  is  given  by 


N2 


(66) 


and  it  is  related  to  the  stability  in  the  circulation  motion  of  the  liquid  (see,  for  example, 
Ref.  17).  Equation  (64)  can  then  be  interpreted  as  providing  a  connection  between  the 
rotation  of  the  Earth  and  the  stability  of  the  induced  motion  within  its  liquid  layer.  For 
this  reason  the  quantity  N,  which  does  not  vanish  everywhere,  should  be  considered  a  more 
appropriate  parameter  than  the  older  Adams-Williamson  condition  in  studying  the  interface 
between  the  liquid  core  and  the  mantle. 


BOUNDARY  CONDITIONS 

We  must  consider  three  types  of  boundary  conditions: 

•  At  the  center  of  the  configuration,  where  r  =  0,  the  solutions  must  be  regular, 
that  is,  representable  by  means  of  infinite  power  series  in  the  radial  distance.  These  series 
can  then  be  used  as  first  approximations  to  initiate  integration.  In  particular,  for  physical 
reasons,  we  must  have  t/(0)  =  F(0)  =  0. 

•  At  any  surface  representing  a  discontinuity  between  two  layers,  the  variables 
U,  R ,  E,  and  H  must  be  continuous  across  them;  the  variables  V  and  F  can  be  discontinuous. 
This  situation  applies  in  particular  to  the  liquid-solid  interface;  however,  between  two  solid 
regions  separated  by  a  discontinuity  we  could  admit  V  to  be  continuous. 
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•  At  the  outer  surface  (r  =  a1 )  we  shall  consider  two  models.  One  model  consists 
of  assuming  a  constant  atmospheric  pressure  (^  1  bar)  at  r  =  aj  to  be  matched  by  the  normal 
stress  E.  The  harmonic  expansion  for  this  type  of  condition  is  as  follows: 

£0(ai)  %  10"3  kbar  =  1  atm, 

En(ai)  =  0,  n  =1,2,  3,  ....  (67) 

)  =  0.  n  =  0,  1,  2,  .... 

A  second  model  for  the  boundary  conditions  is  to  neglect  both  atmospheric  and 
oceanic  layers  and  consider  a  localized  loading  generated  by  the  total  mass  m0  uniformly 
distributed  atop  a  spherical  cap  of  semiaperture  This  generates  a  mass  density 

6(a1(  0o)  =  (1  -  cos  0o) .  (68) 

The  spherical  harmonic  representation  of  this  mass  distribution  is  given  for  n-  0  by 


.1 


5o(qi-  °q)  =  \  I  5(ai,  do>  dv 
-  cos  0o 


mc 


Aira\ 


and  is  given  for  n  >  1  by 


,1 


M<*1>  °o)  =  -  2  I  5(^i ,0o)Pn(v)di> 

' cosd0 


(2  n  +  l)m0 


47raf  (1  -  cos  0Q)  Jcosg0 


f  Pn(v)dv 


(69) 


^0 

47raf  (1  -  cos  60) 


[P n ~ i ( cos  0o)  -  Pn+\ (cos  0o)]  . 


Here  we  have  used  the  facts  that 


(70) 


dP  +1  dP  _  j 

(2n  +  l)Pn(v)  =  ,  n  =  1,2,3,...  (71) 

and 

Pn(l)  =  1,  n  =  0,  1,  2,  ...,  (72) 

as  given  on  pages  17  and  33  of  Ref.  18. 

When  the  spherical  cap  shrinks  to  a  point  (0o  -*  0),  we  get  the  case  of  a  concentrated 
mass  load  m0.  We  can  handle  the  ensuing  indeterminate  form  by  means  of  the  l’Hospital 
rule  and  by  using  the  two  above-mentioned  relationships  to  obtain 
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Pn_  i(v)  P„ +  i(*'>  ..  (dPn*  1  dPn- 1 

v  .\  1  -  v  „_i  \  dv  dv 

=  (2a  +  1)P„(1 )  =  2a  +  1  .  (73) 

Thus 

(2a  +  l)mn 

lim  6n(a^0o)  =  - =  S„(«i)  (74) 

47raj 

is  the  representation  of  the  loading  in  the  limiting  case  (concentrated  loading),  and  this  is 
valid  even  for  a  =  0  (see  Eq.  (69)). 

At  the  surface  of  the  Earth,  the  radial  and  transversal  components  of  the  stress  must 
then  satisfy  the  conditions 


(2n  +  l)m0 

En(al)  =  ~yQ(al  )5„(aj)  =  - - 5 -  70(ai )  -  n  =  0.1,2 . 

4  TOJ 

(75) 

Fn(a\ )  =0,  a  =  0,  1,  2, 


which  have  been  expressed  only  in  the  first-order  approximation  and  with  the  assumption  of 
a  purely  normal  loading. 


Another  boundary  condition  holds  for  the  deformed  potential.  This  is  obtainable  by 
considering,  first,  the  continuity  of  the  deformed  potential  at  the  boundary  (that  is, 

( Rn )  =  Rn ,  where  the  subscript  e  refers  to  the  exterior  space)  and,  second,  the  discontinuity 
suffered  by  the  normal  derivative  of  the  potential  because  of  the  presence  of  a  surface  mass 
distribution  (this  is  the  Gauss  theorem  of  potential  theory): 


dR^ 

bn 


'  -47rG[p0(o1)Un(o1)  +  5„(aj)]  . 


(76) 


In  first-order  approximation,  n  coincides  with  the  radial  direction.  Also,  since  the  exterior 
potential  is  an  infinite  series  in  (l/r)n  +  1,  we  can  write 


Thus,  at  r  =  Oj  we  have 
9  Rr 
bn 


=  <K)  = 


n  +  1 


Rr 
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since  R^)^)  =  R^).  We  get 

,  n  +  1  f  (2  n  +  1  )m0‘ 

R„(a  i)  +~ —  =  4ttG  p0(a1)t/„(ai)  +  - — ^  .  (79) 

1  L  47raf 

In  terms  of  the  normalized  variables,  we  can  write 

H„(« i)  +^rt-IKn(a1)  =  (2n  +  l)-^-  =  <2n  +  l)70(«i>  •  (80) 

1  af 

Thus,  the  conditions  that  hold  at  the  surface  (r  =  a1)  are  expressed  by  Eqs.  (67)  and  (80) 
or  by  Eqs.  (75)  and  (80). 

TOROIDAL  DEFORMATIONS 

When  rotational  symmetry  is  assumed,  as  has  been  the  case  throughout  this  report, 
the  toroidal  displacement  vector  u  has  components  of  the  form 

oo 

u  =  0,  o  =  0,  w  =  Wn(r)dPn/dO .  (81) 

n  =  l 

The  functions  Wn  constitute  a  set  of  unknown  functions  that  will  be  determined  from  the 
equations  of  motion. 

One  can  immediately  verify  that  for  this  type  of  deformation  both  the  dilatation  A 
and  the  perturbed  density  pl  vanish.  As  a  consequence,  the  Poisson  equation  becomes 
V2Vj  =  -4 irGp1  =  0  and  yields 

r2#"  +  2rR'n  -  n(n  +  l)R„  =  0.  (82) 

The  general  solution  is  Rn  =  rn.  However,  a  potential  function  must  vanish  at  infinity, 
which  means  that  we  must  have  Rn  =  0.  This  leads  to  the  vanishing  of  V\  and  of  the 
components  of  its  gradient  gj .  The  Poisson  equation  is  identically  satisfied  and  does  not 
furnish  any  condition. 

By  straightforward  calculations,  it  is  easy  to  realize  that  only  the  longitude-dependent 
component  (the  <*>  component)  of  the  equations  of  motion  gives  a  contribution,  and  it 
amounts  to  the  equation 
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as  a  new  auxiliary  variable.  It  represents  the  radial  factor  of  the  transversal  component  Tr(p 
of  the  stress.  The  normal  form  of  Eqs.  (83)  can  then  be  written  as  the  system 


w'  =  —  W  +  —  Z 
”n  r  ” n  \i  n  ' 


Zi  = 


(n  +2 )(n  -  1)  —  -  p0u>2 


W„ 


—  zn 

r  n 


(85) 


The  boundary  conditions  are 


Z n(ai )  =  0  , 


denoting  the  vanishing  of  the  transversal  stress  at  the  free  surface,  and 


(86) 


Zn(bx)  =  0, 


(87) 


where  bj  is  the  outer  core  radius.  When  p  =  0,  one  can  see  from  Eqs.  (85)  that  Wn  =  0 
and  Zn  =  0. 

Regularity  conditions  at  the  origin  can  be  easily  obtained.  The  series  expansions  for 
n  ¥=  1  are 


-  r2Wn2  + 

~  r^n\  + 

with  arbitrary  Wn2  and 

4 1  =  i(n  +  2)(n-l)pW„2. 

For  n  =  1  one  gets 

W1  =  rWn  +  •••, 

Z i  -  r^Zn  +  '■'» 

where  both  coefficients  are  arbitrary  parameters. 


(88) 


(89) 


(90) 


THE  RESULTS  AND  FUTURE  NUMERICAL  WORK 

We  shall  briefly  summarize  and  discuss  the  results  that  have  been  obtained  in  this 
report.  We  shall  also  express  some  views  on  the  extensive  numerical  work  which  we  plan 
to  undertake  as  a  followup  of  this  theoretical  work. 
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We  start  by  considering  a  rotating  spheroidal  Earth  which  is  in  hydrostatic  equilibrium 
under  the  influence  of  its  self-gravitation  and  its  rotational  acceleration.  We  further  assume 
that  the  rotational  axis  is  fixed  with  respect  to  the  spheroid  and  that  the  angular  velocity  is 
constant.  We  next  consider  linear  perturbing  effects  due  to  the  elastic  property  of  the 
various  layers  comprising  the  interior  of  our  planet.  The  dilatation  of  the  elastic  material 
yielding  to  the  gravity  force  and  centrifugal  force  gives  rise  to  a  perturbed  density  distribu¬ 
tion  within  the  Earth  which  is  obtainable  from  the  continuity  equation  and  which  generates 
a  perturbed  potential  through  the  related  Poisson  equation. 

Elaborating  on  the  equation  of  angular  momentum,  we  ultimately  reach  a  linearized 
version  of  the  Navier-Stokes  equation  and  look  for  solutions  in  the  form  of  spheroidal  and 
toroidal  deformations.  Seismological  evidence  has  corroborated  the  existence  of  such  per¬ 
manent  deformations.  We  have  developed  all  the  variables  in  series  of  Legendre  polynomials, 
and  we  have  ignored  the  terms  due  to  the  Coriolis  acceleration,  since  we  are  assuming  that 
there  is  no  relative  motion  with  respect  to  the  rotating  spheroid.  We  are  limiting  ourselves 
to  studying  the  centrifugal  acceleration  effects.  Further  studies  of  the  oscillatory  motion  of 
the  Earth  will,  however,  require  the  retention  of  the  Coriolis  acceleration. 

The  toroidal  deformations  are  the  easier  of  the  two  cases  to  handle,  because  they  are 
not  associated  with  any  dilation.  Their  study  depends  on  one  second-order  differential 
equation  which  is  related  to  the  longitudinal  component  of  the  Navier-Stokes  equation. 
The  treatment  of  the  spheroidal  deformations  on  the  other  side  gives  rise  to  three  second- 
order  differential  equations:  two  of  them  depend  on  the  radial  and  transversal  components 
of  the  Navier-Stokes  equation;  the  third  one  depends  on  the  Poisson  equation.  The  three 
unknowns  are  the  radial  and  transversal  components  of  the  displacement  vector  field  and 
the  perturbed  potential.  The  reduction  of  the  said  second-order  system  to  its  normal  form, 
a  necessary  feature  when  one  wants  to  perform  any  numerical  evaluation,  introduces  the 
components  of  the  stress  and  the  gravitational  flux  as  auxiliary  variables. 

The  retention  of  the  centrifugal  acceleration  in  the  equations  of  motion  produces 
terms  which  consist  of  the  product  of  Legendre  polynomials  and  of  their  derivatives.  The 
representation  of  the  spheroidal  components  of  such  products  has  constituted  one  of  our 
major  mathematical  problems.  In  Appendix  B  we  elaborate  on  the  well-known  Adams- 
Neumann  product  formula  of  the  Legendre  polynomials  and  succeed  in  expressing  the 
spheroidal  components  of  the  product  of  the  mixed  terms  and  of  the  derivative  terms 
without  introducing  any  new  set  of  harmonic  functions. 

Another  problem  of  an  analytical  character  which  was  encountered  in  the  development 
of  the  equations  of  motion  was  to  provide  the  partial  derivatives  of  the  first  and  second 
order  for  the  unperturbed  potential  with  respect  to  the  radial  distance  and  the  colatitude. 
This  task  was  achieved  quite  simply  in  Appendix  A  by  making  use  of  previously  obtained 
results  pertaining  to  the  equilibrium  configuration  of  a  rotating  planet  which  had  been 
published  in  the  astrophysical  literature.  Our  results  depend  on  the  solution  of  the  Clairaut 
equation  and  have  been  expressed  up  to  the  first  power  of  a  rotational  parameter  which 
characterizes  the  deformation  of  the  original  sphere  into  a  spheroid. 

Of  the  six  differential  equations  comprising  our  normalized  system,  only  two,  the 
equations  expressing  the  derivatives  of  the  stress,  contain  rotational  terms.  Because  of  these 
rotational  terms,  each  equation  representing  the  nth  order  harmonic  of  the  stress  contains 
also  the  harmonics  of  order  n-  2  and  n  +  2  of  some  of  the  other  variables.  Thus  we  have  to 
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deal  with  a  system  simultaneously  relating  three  different  orders  of  harmonics;  this  will 
present  some  difficulty  in  its  numerical  treatment,  since  some  assumption  has  to  be  made 
regarding  the  harmonics  of  order  higher  than  the  order  we  are  solving  or  have  solved  for. 

Our  next  task  will  be  to  numerically  solve  the  system  of  six  differential  equations.  For 
this  purpose  we  plan  to  use  a  Runge-Kutta  integrator  software  package,  since  this  approach 
was  found  satisfactory  by  Longman  19,10]  and  Farrell  ]  1 1  ]  in  dealing  with  the  case  of  a 
nonrotating  Earth.  Most  likely  we  shall  use  DVERK,  which  is  a  fifth-order  integrator 
authored  by  Vemer  [19]  and  which  has  been  discussed  by  Jackson  et  al.  [20] .  We  shall 
be  dealing  with  a  boundary-value  problem  with  conditions  to  be  satisfied  both  at  the  center 
of  the  configuration  and  at  the  free  surface,  and  where  certain  discontinuities  of  the  vari¬ 
ables  can  be  allowed  at  the  interface  between  layers.  Different  Earth  models  will  be  used. 
However,  we  plan  to  devote  more  time  to  the  1066  model  which  was  developed  by  Gilbert 
and  Dziewonski  [21] .  It  consists  of  about  160  data  points  for  the  density  and  the  elastic 
parameters,  and  it  is  considered  to  be  the  most  sophisticated  of  the  existing  models. 

At  the  center  of  the  configuration,  the  displacement  vector  shall  vanish,  and  all  of  the 
six  variables  must  be  regular  functions  of  the  radial  distance.  We  have,  therefore,  obtained 
power-series  solutions  valid  in  the  neighborhood  of  the  origin  and  which  satisfy  the  given 
differential  system  of  equations.  Some  of  the  coefficients  of  the  lowest  powers  appearing 
in  these  series  are  free  parameters  and  must  be  used  in  initiating  the  numerical  solution  and 
in  attempting  to  satisfy  the  other  conditions  at  the  free  surface.  Another  parameter  that 
can  play  a  useful  role  is  the  initial  distance  from  the  origin  where  the  function  evaluations 
are  being  performed  for  the  first  time. 

The  presence  of  the  liquid  outer  core  where  the  rigidity  vanishes  presents  considerable 
complications.  In  evaluating  the  limit  of  the  general  equations  when  the  rigidity  approaches 
zero,  we  Find  that  the  transversal  component  of  the  stress  must  vanish  and,  at  the  same 
time,  that  the  transversal  component  of  the  displacement  can  be  chosen  arbitrarily.  Thus 
these  two  variables  can  suffer  discontinuities  at  the  interface  between  the  inner  and  outer 
cores  and  also  at  the  transition  between  the  outer  core  and  the  mantle.  Further  analysis 
reveals  that  we  are  left  with  four  differential  equations  and  one  algebraic  relationship. 

Our  fundamental  result  in  this  connection  is  that,  due  to  the  adoption  of  a  rotating 
Earth  model,  we  are  not  compelled  to  fall  back  to  the  Adams-Williamson  condition  pertain¬ 
ing  to  the  density  distribution  within  the  outer  core,  which  was  never  completely  accepted 
in  the  past.  We  tend  to  agree  with  other  researchers  to  pom:  out  that  a  better  parameter  to 
be  used  at  this  interface  is  the  Brunt-Vaisala  frequency  for  the  circulation  motion  of  the 
fluid.  Another  approach  which  we  shall  also  consider,  and  which  might  have  more  realistic 
physical  implications ,  is  the  adoption  of  Moloc*  mskii’s  analytic  theory  [6]  to  represent  the 
motion  within  the  liquid  layer. 

At  the  other  surfaces  of  discontinuity  within  the  mantle,  we  shall  suppose  that  all  six 
variables  are  continuous  unless  numerical  considerations  will  otherwise  impose  upon  us  a 
different  outlook.  At  the  free  surface,  we  have  three  conditions.  One  is  obtained  by 
applying  Gauss  theorem  to  the  discontinuity  suffered  by  the  mass  distribution.  The  second 
condition  is  the  vanishing  of  the  tangential  stress.  The  third  condition  is  the  matching  of 
the  radial  stress  with  the  loading  conditions.  We  are  assuming  a  concentrated  load  and  have 
followed  Longman’s  approach  in  obtaining  its  analytic  representation  as  a  limit  of  a  load 
uniformly  distributed  upon  a  small  spheroidal  cap  when  the  aperture  of  this  cap  goes 
to  zero. 
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Starting  from  the  fundamental  system  of  equations  corresponding  to  the  zero-order 
harmonic,  we  can  recursively  obtain  the  Love  numbers  and  the  load  numbers  of  any  order. 

It  will  be  interesting  at  this  stage  to  study  the  variation  of  these  two  sets  of  numbers  with 
the  particular  Earth  model  used.  Speculation  by  both  Farrell  [11]  and  Varga  [22]  has 
been  that  the  Love  numbers  of  large  orders  are  not  much  influenced  by  the  structure  of  the 
core  but  depend  essentially  on  the  values  that  the  variables  assume  at  the  base  of  the  mantle. 

Once  the  load  numbers  have  been  ascertained  up  to  a  certain  finite  order,  one  has  to 
consider  the  sum  of  the  infinite  Legendre  polynomial  expansions  to  obtain  the  spheroidal 
deformations  and  other  pertinent  information  about  the  tilt  and  tidal  gravity.  For  this 
purpose,  we  shall  elaborate  upon  certain  procedures  followed  by  Farrell  in  Ref.  11  that  deal 
with  the  nonrotating  Earth.  These  procedures  entail:  the  use  of  certain  algorithms  to  im¬ 
prove  the  rapidity  of  convergence  of  the  series  involved,  the  development  of  an  asymptotic 
theory  to  ascertain  the  appropriate  location  where  the  truncation  of  the  series  should  take 
place,  the  sum  of  the  Poisson  series  of  Legendre  polynomials  for  values  of  the  argument  less 
than  1  to  provide  the  value  of  the  deformation  at  a  certain  angular  distance  from  the  loca¬ 
tion  of  the  applied  load,  and  the  use  of  the  well-known  Boussinesq  [23]  solution  for  an 
elastic  flat  plate  to  approximate  the  values  that  the  variables  assume  at  points  in  the 
immediate  neighborhood  of  the  applied  load. 
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Appendix  A 

DERIVATIVES  OF  THE  UNPERTURBED  POTENTIAL 


We  must  obtain  the  partial  derivatives  of  the  unperturbed  potential  V0(r,  9)  pertaining 
to  hydrostatic  equilibrium  with  respect  to  both  variables  r  r  .id  0.  For  this  purpose,  we 
refer  to  formula  (27)  in  Ref.  Al,  in  which  reference  the  said  potential  was  determined  in 
terms  of  the  mean  radius  a  which  characterizes  an  equipotential  surface.  More  specifically,  let 

r  =  a[l  +  q/‘21(a)P2(cos  0)]  +  (Al) 

be  the  equation  of  an  equipotential  surface  to  first-order  terms  in  the  rotational  parameter 

9  =  =  : ~  r  : »  (A2> 

iGmj  4ffGp0(a1) 

where  is  the  mean  radius  of  the  outermost  surface  and  where  the  nondimensional 
function  f2i,  a  solution  of  the  Clairaut  equation,  represents  the  spheroidal  deformation.  In 
Ref.  Al  the  potential  V0(r,  0)  was  determined  along  an  equipotential  surface  r(a,  9)  as  a 
function  of  a : 


V0[r(a,  0);  0]  H  »Ma)  ■  (A3) 

We  use  this  function  0o(a)  to  calculate  the  required  derivatives.  This  can  be  achieved  most 
expeditiously  by  recalling  a  well-known  expansion  theorem  due  to  Lagrange  (see,  for 
example.  Ref.  A2),  which  states:  if  the  variables  a  and  r  are  implicitly  related  according 
to  the  relation 


then 


a  =  r  +  q<t>(a) , 


dF(r) 

F(a)  =  F(r)  +  q<p(r)  +  •••, 

where  F  stands  for  a  generic  function.  Equation  (Al)  reveals  that 

0(a)  =  -af21(a)P2(cos  9) , 

and  this  leads  to 


tf'ofa)  =  +  9 [_r/2l(r)^>2(cos  0)1 


4jtG 


900 

dr 


(A4) 


(A5) 


|  r2p0(r)  +  f  p0(r)rdr 
Jr 

+  ^-  u2r2A*(r)P2{cos  9)  +  cu2r2  +  •••. 


(A6) 
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Here  we  have  set 


A*(r) 


Pp(r) 

Po(ri) 


hi(r) , 


a  nondimensional  function  of  r,  and  have  used  the  formula 


Po 


=  7-  (p0  -  Po)  . 


(A7) 


(A8) 


where  primes  as  usual  denote  derivatives  with  respect  to  r.  By  differentiation  of  Eq.  (A6) 
with  respect  to  r  and  6 ,  we  get 


3  Vn  a  o  „  1 

—  =  --o  irGrp0(r)  +  —  w2r  +  u>2rB*(r)P2(cos  9) , 


dr 


dV, 


-  =  —  cj2r2A*(r) 

30  3  u  r  A  (r)  dQ  , 


d2V, 


dr2 


°  =  4  irG(2p0  -  3p0)  +  4  co2  -4-  co2C*(r)P2(cos  0) , 


drdO  3  W  B  (  }  30  ’ 


92V0  =  1 

302  3 


o  o  *  d2p2 

w2r2A*(r)  — — 
302 


(A9) 


The  following  additional  notation  has  been  used 


B*(r)  =  (3D0  +  r)2l  ~  D/kW, 

Po(rl) 

C*(r)  =  Ld0  +  2t721  -  8  -  3  f21(r) . 

PoOt)  \  Po/ 


(A10) 


These  are  nondimensional  expressions  which  depend  on  the  solution  of  the  Clairaut 
equation.  We  have  also  set 


D0  =  Pq/Pq  ,  ??  =  rf'/f 

and  have  eliminated  the  terms  containing  rj'  by  means  of  the  Radau  equation 

rr?'  =  6  +  ??(1  -  t?)  -  6D0(1  +  rj) 


(All) 


(A12) 
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(see  Ref.  A3).  It  is  easy  to  verify  that 

B*  =  2A*  +  r(A*)' , 
-C*  =  B*  +  r(B*)' . 


(A13) 


These  relations  can  be  used  to  generate  B*  and  C*  if  one  wants  to  avoid  numerically  differ¬ 
entiating  p0.  For  computational  purposes,  it  is  appropriate  to  note  that  the  limit  of  the 
three  functions  A*(r),  B*(r),  and  C*(r)  when  r  approaches  zero  is  also  zero. 
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Appendix  B 


SPHEROIDAL  COMPONENTS  FOR  THE  PRODUCT  OF 
TWO  LEGENDRE  POLYNOMIALS  AND  THEIR  DERIVATIVES 


We  start  from  the  well-known  Adams-Neumann  formula  which  was  proved  independ¬ 
ently  by  each  of  these  authors  in  1878  [B1,B2]  (see  Refs.  B3,  B4,  and  18).  The  said 
formula  expresses  the  product  of  two  Legendre  polynomials  as  a  sum  of  the  same 
polynomials: 


\  — 


<7 

(cos  (9)?,  (cos  0)  =  J2  B(p,q-,r)Pp+q.2r  (cos  6). 

r= o 


(Bl) 


Here,  p  >  q  are  positive  integers,  and  the  B  coefficients  are  given  by 

rwn  „  .  _  Ap-rArAq-r  2p  +  2q  -  4r  +  1 

B(P,  q;  r )  = - 

with 


Ap+q-r  2p  +  2q  -  2r  +  1  ’ 


1  •  3- 5  ...  (2r-  1) 

Ar  - - ~i - »  r  =  1,  2,  ..., 


Ao  ~  1  • 


(B2) 


(B3) 


Note  that  p  and  q  appear  symmetrically  within  B. 

We  differentiate  Eq.  (Bl)  twice  with  respect  to  9  and  eliminate  the  second  derivatives 
by  using  the  well-known  result  (see  Eq.  (28)) 


_  =  -(cot*)  —  -  p(p  +  l)Pp. 


(B4) 


After  some  reordering  of  the  terms,  we  get 


Wp  Wq 

2wW~  (cot0) 


30  (PpP^ 


dP. 


r=0 


p*q~2r 

~de 


=  lp(p  +  1)  +  q(q  +  l)]PpPq 


q 

(P  +  q~  2r)(p  +  q-2r  +  l)BPp+9_2r.  (B5) 

r  =  0 

32 
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The  bracketed  term  which  appears  in  the  left-hand  side  vanishes,  because  it  is  a  differential 
consequence  of  Eq.  (Bl).  We  next  proceed  to  transform  the  product  PpPq  in  the  right-hand 
side  by  using  again  Eq.  (Bl )  and  get  f  ' 


2  dfT  JF  =  2J  tp(p  +  D  +  Q(Q  +  1)  -  (P  +  Q  -  2 r)(p  +  q  ~  2r  +  1  )]BPp+q.2r.  (B6) 


We  can  then  write 


where 


9pP  3 Pq  « 

30  30  =  C(P’  9',  r)Pp+9-2r> 

r  =  0 


C(p,  q;  r)  =  [r(2p  +  2q  -  2r  +  1)  -  pq]B(p,  q;  r) . 


Note  that  the  coefficient  C  is  symmetric  with  respect  to  p  and  q.  Let  us  now  tum  to  the 
product  P  bP/bO  and  assume  that  it  can  be  represented  as 

55-  P,  -  ^  fl(p,  r>  -^5— .  CB! 

r=0 

It  is  a  matter  then  of  determining  an  algebraic  expression  for  D.  We  differentiate  Eq.  (B9) 
once  with  respect  to  9  and  eliminate  the  second  derivatives  through  Eq.  (B4)  and  write 


■ 

(P  +  q-  2r)(p  +  q-  2r+l)DPp+(. 


bPp  bPq  bPp  Q  dPp  +  q-2r 

=  P(P  +  l)PpPq  -  gg  W  *  (COt  0)  W  Pq  ~  I  D  — 

L  r  =  0 


.  (BIO) 


The  term  within  brackets  vanishes  because  of  Eq.  (B9).  We  next  replace  both  products  in 
the  right-hand  side  by  means  of  Eqs.  (Bl)  and  (B7).  We  equate  the  coefficients  of  the 
polynomials  of  like  order  and  get 


(p  +  q  -  2r)(p  +  q  -  2r  +  1)D  =  p(p  +  l)B  -  C. 


(Bll) 


Recalling  the  expression  for  C  given  by  Eq.  (B8),  we  can  finally  write 


„  P(P  +  1)  +  pq  -  r(2p  +  2q  -  2r  + 1) 

D(P’  r)  = - (P  +"q  ~  2r)(p  +  ,'r2rTi) -  B(p,  q;  r) . 


(B12) 
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This  expression  is  not  symmetric  in  p  and  q .  Naturally 


Also,  since 


we  must  have 


bPq  y  i  dPp+q-2r 

PP  ~b6  =  L  Dto-P'-r)  - 


K  bpc  a 

— -  p  +  p  — -  =  _  ip  p  i 

be  *  p  be  be  K  p  ’ 


D(p,q\r)  +  D(q,  p;  r)  =  B(p,q\r), 


(B13) 


(B14) 


(B15) 


and  this  can  be  directly  verified.  When  q  =  0,  the  summations  appearing  in  Eqs.  (Bl),  (B7), 
and  (B9)  reduce  to  only  one  term,  which  is  due  to  r  =  0.  The  limiting  values  for  the 
coefficients  are 


(B16) 


B(p,  0;  0)  =  1 ,  C(p,  0;  0)  =  0  , 

D(p,  0;0)  =  1,  D(0,  p;  0)  =  0. 

We  have  used  the  previously  obtained  formulas  to  evaluate  the  products 

95,  85,  85,  Sffc 

”  2 '  88  88  ’  88  2 '  "  88 


which  appear  in  the  fundamental  equations  of  motion  for  a  generic  value  of  n.  The  results 
are  presented  in  Table  Bl.  The  terms  multiplied  by  the  factor  (1  -  5ln),  where  6  * .  the 
Kronecker  delta,  will  not  appear  when  n  =  1. 


Table  Bl  —  Resolution  of  Products  of  Legendre  Polynomials 
Into  Spheroidal  Components 


3 (n  +  l)(n  +  2) 

2  n(n  +  1) 

3 (n  -  1  )n 

~*n*2 

(2n  +  l)(2n  +  3)  ‘n  +  2 

(2n  -  l)(2n  +  3)  * 

"  '  (2n  -  l)(2n  +  1)  Pn~ 2 

1  *P2 

n{n  +  1)  («  +  2) 

n(n  +  1) 

n(n-l)(n  +  l) 

3  be  be 

(2n  +  l)(2n  +  3)  />"  +  2 

(2n  -  l)(2n  +  3) 

*"  (2n  -  l)(2n  +  1)  '"-2 

bPn 

3n(n  +  1)  bPn+2 

2(n2  +  n  -  3) 

9Pn  |  3n(n  +  1)(1  -  6ln)  bPn_2 

2  be  p 2 

(2 n  *  l)(2n  +  3)  b6 

(2n  -  1) (2n  +  3) 

be  '  (2n  -  l)(2n  +  1)  be 

ip  ^ 

n  +  1  **Pn  +  2 

,  1 

bPn  n(l-6ln)  bPn_2 

3  "  be 

(2n  +  l)(2n  +  3)  be 

(2 n  -  l)(2n  +  3) 

be  (2n  -  1)(2tj  +  1)  be 
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